{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "# Gaussian Mixture Model Selection\n",
        "\n",
        "This example shows that model selection can be performed with\n",
        "Gaussian Mixture Models using information-theoretic criteria (BIC).\n",
        "Model selection concerns both the covariance type\n",
        "and the number of components in the model.\n",
        "In that case, AIC also provides the right result (not shown to save time),\n",
        "but BIC is better suited if the problem is to identify the right model.\n",
        "Unlike Bayesian procedures, such inferences are prior-free.\n",
        "\n",
        "In that case, the model with 2 components and full covariance\n",
        "(which corresponds to the true generative model) is selected.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {},
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import itertools\n",
        "\n",
        "from scipy import linalg\n",
        "import matplotlib.pyplot as plt\n",
        "import matplotlib as mpl\n",
        "\n",
        "from sklearn import mixture"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Number of samples per component\n",
        "n_samples = 500\n",
        "\n",
        "# Generate random sample, two components\n",
        "np.random.seed(0)\n",
        "X = np.r_[\n",
        "    np.random.randn(n_samples, 1),\n",
        "    0.7 * np.random.randn(n_samples, 1) + 3,\n",
        "]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<AxesSubplot:ylabel='Count'>"
            ]
          },
          "execution_count": 27,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD4CAYAAAAD6PrjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAATXUlEQVR4nO3df5Af9X3f8ecLZMA/axgORj6dKtLK1JhxbM+ZYmgztjE1TTzIyYQgJnbUhkRui107SW1DmYbpdDTD1JnUnaZOrNoU0lAwcXAgdoJRZGwmYxt82NhGCIIm1OgsBZ3DZOzWHajg3T9utf5y3Emnk/a7X+n7fMzcfL/72f3uvU5I92J3v/v5pqqQJAnghL4DSJJGh6UgSWpZCpKklqUgSWpZCpKk1qq+AxyJ008/vdatW9d3DEk6pjzwwAPfr6qJxdYd06Wwbt06ZmZm+o4hSceUJN9dap2njyRJLUtBktSyFCRJLUtBktSyFCRJLUtBktSyFCRJLUtBktSyFCRJLUtB0hGZnFpLks6+JqfW9v0jjpVjepoLSf3bM7ubyz/xlc72/+n3XtDZvvVCHilIklqWgiSpZSlIklqdlUKSG5LsS/LQgvH3J3k0yY4k/3Fg/Joku5p17+gqlyRpaV1eaL4R+B3g9w8MJHkrsAF4XVU9neSMZvwcYCPwWuBVwJ8neXVVPdthPknSAp0dKVTVvcBTC4b/JXB9VT3dbLOvGd8A3FpVT1fV48Au4LyuskmSFjfsawqvBv5xkvuSfDnJm5rxSWD3wHazzdgLJNmcZCbJzNzcXMdxJWm8DLsUVgGnAucDHwJuSxIgi2xbi+2gqrZW1XRVTU9MLPoRo5KkFRp2KcwCt9e8+4HngNOb8amB7dYAe4acTZLG3rBL4Y+BtwEkeTVwEvB94E5gY5KTk5wFrAfuH3I2SRp7nb37KMktwFuA05PMAtcBNwA3NG9TfQbYVFUF7EhyG/AwsB+4ynceSdLwdVYKVXXFEqvevcT2W4AtXeWRJB2adzRLklqWgiSpZSlIklqWgiSpZSlIklqWgiSpZSlIklqWgiSpZSlIklqWgiSpZSlIklqWgiSpZSlIklqWgiSpZSlIklqWgiSp1VkpJLkhyb7mU9YWrvs3SSrJ6QNj1yTZleTRJO/oKpckaWldHincCFyycDDJFHAx8MTA2DnARuC1zWs+nuTEDrNJkhbRWSlU1b3AU4us+k/Ah4EaGNsA3FpVT1fV48Au4LyuskmSFjfUawpJLgW+V1XfWrBqEtg9sDzbjC22j81JZpLMzM3NdZRUksbT0EohyUuAa4HfXGz1ImO1yBhVtbWqpqtqemJi4mhGlKSxt2qI3+vvAWcB30oCsAb4RpLzmD8ymBrYdg2wZ4jZJEkM8Uihqr5TVWdU1bqqWsd8Ebyxqv4auBPYmOTkJGcB64H7h5VNkjSvy7ek3gJ8FTg7yWySK5fatqp2ALcBDwN3AVdV1bNdZZMkLa6z00dVdcUh1q9bsLwF2NJVHknSoXlHsySpZSlIklqWgiSpZSlIklqWgiSpZSlIklqWgiSpZSlIklqWgqTRdsIqknTyNTm1tu+fbuQMc0I8STp8z+3n8k98pZNdf/q9F3Sy32OZRwqSpJalIElqWQqSpJalIElqWQqSpJalIElqdfnJazck2ZfkoYGxjyZ5JMm3k3w2ySsH1l2TZFeSR5O8o6tckqSldXmkcCNwyYKxbcC5VfU64C+BawCSnANsBF7bvObjSU7sMJskaRGdlUJV3Qs8tWDs7qra3yx+DVjTPN8A3FpVT1fV48Au4LyuskmSFtfnNYVfBv6seT4J7B5YN9uMvUCSzUlmkszMzc11HFGSxksvpZDkWmA/cPOBoUU2q8VeW1Vbq2q6qqYnJia6iihJY2nocx8l2QS8E7ioqg784p8FpgY2WwPsGXY2SRp3Qz1SSHIJ8BHg0qr60cCqO4GNSU5OchawHrh/mNkkSR0eKSS5BXgLcHqSWeA65t9tdDKwLQnA16rqX1TVjiS3AQ8zf1rpqqp6tqtskqTFdVYKVXXFIsOfOsj2W4AtXeXRaJicWsue2d2H3nCFXrVmiu/tfqKz/UvHOz9PQUO1Z3Z3Z3Pjg/PjL6XrMtbxw1KQxkCXZWwRH1+c+0iS1LIUJEktS0GS1LIUJEktS0GS1LIUJEktS0GS1LIUJEktS0GS1LIUJEktS0GS1LIUJEktS0GS1LIUJEmtzkohyQ1J9iV5aGDstCTbkjzWPJ46sO6aJLuSPJrkHV3lkiQtrcsjhRuBSxaMXQ1sr6r1wPZmmSTnABuB1zav+XiSEzvMJklaRGelUFX3Ak8tGN4A3NQ8vwl418D4rVX1dFU9DuwCzusqmyRpccO+pnBmVe0FaB7PaMYngcHPCpxtxl4gyeYkM0lm5ubmOg0rSeNmVC40Z5GxWmzDqtpaVdNVNT0xMdFxLB1zTlhFkk6+JqfW9v3TSZ0b9mc0P5lkdVXtTbIa2NeMzwJTA9utAfYMOZuOB8/t97OIpSOwrCOFJBcuZ2wZ7gQ2Nc83AXcMjG9McnKSs4D1wP0r2L8k6Qgs9/TRf1nmWCvJLcBXgbOTzCa5ErgeuDjJY8DFzTJVtQO4DXgYuAu4qqqeXWY2HWWTU2s7OwUjabQd9PRRkjcDFwATSX59YNUrgIO+ZbSqrlhi1UVLbL8F2HKwfWo49szu9hSMNKYOdU3hJOBlzXYvHxj/AfDzXYWSJPXjoKVQVV8Gvpzkxqr67pAySZJ6stx3H52cZCuwbvA1VfW2LkJJkvqx3FL4Q+D3gE8CXgCWpOPUckthf1X9bqdJJEm9W+5bUv8kyb9KsrqZ6fS0JKd1mkySNHTLPVI4cMPZhwbGCviJoxtHktSnZZVCVZ3VdRBJUv+WVQpJfmmx8ar6/aMbR5LUp+WePnrTwPNTmL8r+RuApSBJx5Hlnj56/+Bykr8D/I9OEkmSerPSz1P4EfMzmUqSjiPLvabwJ/z4Q29OBF7D/KymkqTjyHKvKfzWwPP9wHeraraDPJKkHi3r9FEzMd4jzM+UeirwTJehJEn9WO4nr/0C85+EdhnwC8B9SVY8dXaSX0uyI8lDSW5Jckpzl/S2JI81j6eudP+SpJVZ7oXma4E3VdWmqvol4Dzg363kGyaZBP41MF1V5zJ/jWIjcDWwvarWA9ubZUnSEC23FE6oqn0Dy39zGK9dzCrgxUlWAS8B9gAbgJua9TcB7zqC/UuSVmC5F5rvSvIF4JZm+XLgT1fyDavqe0l+C3gC+L/A3VV1d5Izq2pvs83eJGesZP+SpJU71Gc0/33gzKr6UJKfA/4REOCrwM0r+YbNtYINwFnA3wJ/mOTdh/H6zcBmgLVr164kgiRpCYc6BfQx4IcAVXV7Vf16Vf0a80cJH1vh93w78HhVzVXV/wNuBy4AnkyyGqB53LfYi6tqa1VNV9X0xMTECiNIkhZzqFJYV1XfXjhYVTPMfzTnSjwBnJ/kJUnC/DxKO4E7+fEU3ZuAO1a4f0nSCh3qmsIpB1n34pV8w6q6L8lnmJ9Qbz/wTWAr8DLgtiRXMl8cl61k/5KklTtUKXw9ya9W1X8bHGx+cT+w0m9aVdcB1y0Yfpr5owZJUk8OVQofBD6b5Bf5cQlMAycBP9thLklSDw5aClX1JHBBkrcC5zbDn6+qL3aeTJI0dMv9PIV7gHs6ziJJ6tmR3JUsSTrOWAqSpJalIElqWQqSpJalIElqWQrScp2wiiSdfE1OObmjRsNyp86W9Nx+Lv/EVzrZ9affe0En+5UOl0cKkqSWpSBJalkKkqSWpSBJalkKkqSWpSBJavVSCklemeQzSR5JsjPJm5OclmRbkseax1P7yCZJ46yvI4X/DNxVVf8A+EnmP6P5amB7Va0HtjfL0njo8Ma4+Y9Cl5Zn6DevJXkF8FPAPwOoqmeAZ5JsAN7SbHYT8CXgI8POJ/WiwxvjwJvjtHx9HCn8BDAH/Pck30zyySQvBc6sqr0AzeMZi704yeYkM0lm5ubmhpdaksZAH6WwCngj8LtV9Qbg/3AYp4qqamtVTVfV9MTERFcZJWks9VEKs8BsVd3XLH+G+ZJ4MslqgOZxXw/ZJGmsDb0Uquqvgd1Jzm6GLgIeBu4ENjVjm4A7hp1NksZdX7Okvh+4OclJwF8B/5z5grotyZXAE8BlPWWTpLHVSylU1YPA9CKrLhpyFEnSAO9oliS1LIVj0OTUWm90ktQJP3ntGLRndrc3OknqhEcKkqSWpSBJalkKkqSWpSBJalkKksZXh1OWT06t7funWxHffSRpfHU4Zfmx+i4+jxQkSS1LQZLUshQkSS1LQZLUshQkSS1LQZLUshQkSa3eSiHJiUm+meRzzfJpSbYleax5PLWvbJI0rvo8UvgAsHNg+Wpge1WtB7Y3y5KkIeqlFJKsAX4G+OTA8Abgpub5TcC7hhxLksZeX0cKHwM+DDw3MHZmVe0FaB7PWOyFSTYnmUkyMzc313lQSRonQy+FJO8E9lXVAyt5fVVtrarpqpqemJg4yukkabz1MSHehcClSX4aOAV4RZI/AJ5Msrqq9iZZDezrIZskjbWhHylU1TVVtaaq1gEbgS9W1buBO4FNzWabgDuGnU2Sxt0o3adwPXBxkseAi5tlSdIQ9fp5ClX1JeBLzfO/AS7qM48kjbtROlKQJPXMUpAktSwFSVLLUpAktSwFSVLLUpAktSwFSVLLUpAktSwFSVLLUpAktSwFSVLLUpAktSwFSVLLUpAktSwFSVLLUpAktYZeCkmmktyTZGeSHUk+0IyflmRbkseax1OHnU2Sxl0fRwr7gd+oqtcA5wNXJTkHuBrYXlXrge3NsiRpiIZeClW1t6q+0Tz/IbATmAQ2ADc1m90EvGvY2SRp3PV6TSHJOuANwH3AmVW1F+aLAzhjiddsTjKTZGZubm5oWSVpHPRWCkleBvwR8MGq+sFyX1dVW6tquqqmJyYmugt4hCan1pKkky9J6sqqPr5pkhcxXwg3V9XtzfCTSVZX1d4kq4F9fWQ7WvbM7ubyT3ylk31/+r0XdLJfSerj3UcBPgXsrKrfHlh1J7Cpeb4JuGPY2SRp3PVxpHAh8B7gO0kebMb+LXA9cFuSK4EngMt6yCZJY23opVBVfwEsdWL8omFmkSQ9n3c0S5JaloIkqWUpSJJaloIkdeGEVZ3dq5SEyam1ncTu5T4FSTruPbe/s3uVoLv7lTxSkCS1LAVJUstSkCS1xroUnLROkp5vrC80O2mdJD3fWB8pSJKez1KQJLUsBUlSy1KQJLUsBUlSy1KQJLVGrhSSXJLk0SS7klzddx5JGicjVQpJTgT+K/BPgXOAK5Kc028qSRofI1UKwHnArqr6q6p6BrgV2NBzJkkaG6mqvjO0kvw8cElV/Uqz/B7gH1bV+wa22QxsbhbPBh5dYnenA9/vMO6RGOVsMNr5RjkbjHa+Uc4Go51vlLPB4ef7u1U1sdiKUZvmYrFJg57XWlW1Fdh6yB0lM1U1fbSCHU2jnA1GO98oZ4PRzjfK2WC0841yNji6+Ubt9NEsMDWwvAbY01MWSRo7o1YKXwfWJzkryUnARuDOnjNJ0tgYqdNHVbU/yfuALwAnAjdU1Y4V7u6Qp5h6NMrZYLTzjXI2GO18o5wNRjvfKGeDo5hvpC40S5L6NWqnjyRJPbIUJEmt47YUkvyHJN9O8mCSu5O8qu9Mg5J8NMkjTcbPJnll35kOSHJZkh1JnksyMm/DG+UpUJLckGRfkof6zrJQkqkk9yTZ2fx3/UDfmQ5IckqS+5N8q8n27/vOtFCSE5N8M8nn+s6yUJL/leQ7ze+5maOxz+O2FICPVtXrqur1wOeA3+w5z0LbgHOr6nXAXwLX9Jxn0EPAzwH39h3kgGNgCpQbgUv6DrGE/cBvVNVrgPOBq0boz+5p4G1V9ZPA64FLkpzfb6QX+ACws+8QB/HWqnr98XqfwlFTVT8YWHwpC26C61tV3V1V+5vFrzF/T8ZIqKqdVbXUneJ9GekpUKrqXuCpvnMspqr2VtU3muc/ZP4X3GS/qebVvP/dLL6o+RqZf6tJ1gA/A3yy7yzDctyWAkCSLUl2A7/I6B0pDPpl4M/6DjHiJoHdA8uzjMgvtmNJknXAG4D7eo7Sak7PPAjsA7ZV1chkAz4GfBh4ruccSyng7iQPNFMAHbFjuhSS/HmShxb52gBQVddW1RRwM/C+g+9t+Pmaba5l/vD+5lHLNmIOOQWKDi7Jy4A/Aj644Ei6V1X1bHOadw1wXpJze44EQJJ3Avuq6oG+sxzEhVX1RuZPq16V5KeOdIcjdfPa4aqqty9z0/8JfB64rsM4L3CofEk2Ae8ELqoh3zByGH92o8IpUI5AkhcxXwg3V9XtfedZTFX9bZIvMX9tZhQu2F8IXJrkp4FTgFck+YOqenfPuVpVtad53Jfks8yfZj2ia4HH9JHCwSRZP7B4KfBIX1kWk+QS4CPApVX1o77zHAOcAmWFkgT4FLCzqn677zyDkkwceOddkhcDb2dE/q1W1TVVtaaq1jH/9+2Lo1QISV6a5OUHngP/hKNQpsdtKQDXN6dDvs38H9bIvA2v8TvAy4FtzdvJfq/vQAck+dkks8Cbgc8n+ULfmZqL8gemQNkJ3HYEU6AcdUluAb4KnJ1kNsmVfWcacCHwHuBtzd+1B5v/+x0Fq4F7mn+nX2f+msLIvfVzRJ0J/EWSbwH3A5+vqruOdKdOcyFJah3PRwqSpMNkKUiSWpaCJKllKUiSWpaCJKllKUiSWpaCJKn1/wErwFBo3T38kAAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "import seaborn as sns\n",
        "\n",
        "sns.histplot(x=X[:, 0])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "metadata": {},
      "outputs": [],
      "source": [
        "lowest_bic = np.infty\n",
        "bic = []\n",
        "n_components_range = range(1, 7)\n",
        "cv_types = [\"spherical\", \"tied\", \"diag\", \"full\"]\n",
        "for cv_type in cv_types:\n",
        "    for n_components in n_components_range:\n",
        "        # Fit a Gaussian mixture with EM\n",
        "        gmm = mixture.GaussianMixture(\n",
        "            n_components=n_components, covariance_type=cv_type\n",
        "        )\n",
        "        gmm.fit(X)\n",
        "        bic.append(gmm.bic(X))\n",
        "        if bic[-1] < lowest_bic:\n",
        "            lowest_bic = bic[-1]\n",
        "            best_gmm = gmm"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x7f9df32f98e0>"
            ]
          },
          "execution_count": 29,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAADRCAYAAADhe+qZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhyUlEQVR4nO3de5xXVb3/8dcbJFG5KXD6yXX0JCaichNREwn5iSdNrTS01Mxf4o1KTY+anRN6xCytX2qJmRiaeSG8JYJiB0g9idyvKYI6/hzlIOJBwQQFP78/9oK+wHeYma8D39nD+/l4fB+zZ+211v7s7eUza6/13VsRgZmZmTV8TcodgJmZmdWOk7aZmVlOOGmbmZnlhJO2mZlZTjhpm5mZ5YSTtpmZWU44aZtZoyNpjKTralm3UtLg7R2TWX1w0jarpfQ/9w8lrZH0P5KekNS5YP9miULSZySNkLRE0gep/V2SKspyAmaWe07aZnXz5YhoAewNLAdu3UbdccCJwDeA1sAhwCzgmO0d5JaUKct/7+U8tllj4/+QzEoQEWvJknL3YvvT7db/DZwUETMiYn1EvBcRv46I0dW0uULSm5JWS1os6ZhU3lTSDyW9kvbN2jjCl3SEpBmS3ks/jyjob6qkkZL+C/g7sK+kz0t6WtK76Rhfr+4cU/ufSJqe+n9M0l4F+/tL+qukVZLmSRq4rWMX6b9S0uWS5qc7EaMlfVbSxHSef5a0Z0H9EyUtSsebKumAgn29JM1O7R4Emm9xrBMkzU1t/yrp4OrO26whc9I2K4Gk3YGhwLRqqgwGpkfEG7Xsb39gOHBoRLQEhgCVafelwOnAl4BWwDnA31MCfQK4BWgL/AJ4QlLbgq7PBIYBLYEVwNPAfcA/pT5vk3TgNkI7Kx2vA7A+HQtJHdOxrwP2Ai4DHpLUvppjv15N/18j++OmG/BlYCLwQ6Ad2f+fvpeO1w24H7gYaA9MAB5PUxCfAR4Ffp9i+WPql9S2N3AXcF66Tr8B/iRp122ct1mD5KRtVjePSloFvE+WbG6spl5bYFkd+t0A7Ap0l9QsIioj4pW07zvAjyJicWTmRcRK4HhgSUT8Po3k7wdeIkt+G42JiEURsR44DqiMiN+l+rOBh4BTthHX7yNiYUR8APwb8HVJTYEzgAkRMSEiPomIp4GZZH9YbHXsiPi4mv5vjYjlEfEm8CzwQkTMiYh1wCNAr1RvKPBERDyd+roJ2A04AugPNAN+GREfR8Q4YEbBMc4FfhMRL0TEhoi4G1iX2pnlipO2Wd2cHBFtyBLscOAvkv5XkXoryea9ayUilpKNIkcAb0t6QFKHtLsz8EqRZh3YegT7OtCx4PfCkX5X4LB0i3hV+uPjm0Cx+Iu1f50sObZLfZ26RV9fYPNzrs1dhuUF2x8W+b1F2t7sXCPik9R/x7Tvzdj87UeF16Ur8IMtYu2c2pnlipO2WQnSiO1hshHyF4pU+TPQT1KnOvR5X0R8gSzJBPDTtOsN4J+LNHkr1S3UBXizsNuC7TeAv0REm4JPi4i4YBthdS7Y7gJ8DLyT+vr9Fn3tERE3VHPsT2uzc5WkFNubZHc0Oqaywlg3egMYuUWsu6c7E2a54qRtVoK0IvokYE/gxS33R8SfyeaPH5HUR9IuklpKOl/SOUX621/SoDTPupZslLkh7b4T+A9J+6XjHpzmrScA3SR9I/U/lGxh3Phqwh6f6p8pqVn6HFq4oKuIMyR1T3P41wLjImIDcC/wZUlD0kK55pIG1uWPlDoaCxwv6RhJzYAfkN3i/ivwPNl8+/fSdfgq0K+g7W+B8yUdlq7fHpKOl9RyO8Vqtt04aZvVzeOS1pDNaY8EvhURi6qpewpZYn0QeA9YCPQlG4VvaVfgBrJR7H+TLRT7Ydr3C7KkNSkddzSwW5rXPoEsga0E/hU4ISLeKRZMRKwGjgVOIxu5/jfZaH5bC7J+D4xJdZuTFoalBXYnpRhXkI1mL2c7/T8lIhaTzaPfSnaNvkz29buPIuIj4KvA2cD/kM1/P1zQdibZvPav0v6lqa5Z7mjzaSAzs4ykqcC9EXFnuWMxs4xH2mZmZjnhpG1mZpYTvj1uZmaWEx5pm5mZ5YSTtpmZWU7sUu4AatKuXbuoqKiot/5mzXqrpHYtDtqz5kpFNH+vWUntun4yr6R2fLZPae3MzKxBmDVr1jsR0b7YvgaftCsqKpg5c2a99SddU1K7vo+fXFK7zz3RueZKRfz2w7Y1VyrmB/V3rczMbMeTVN0Ldnx73MzMLC+ctM3MzHLCSdvMzCwnGvyctpmZNUwff/wxVVVVrF27ttyh5FLz5s3p1KkTzZrVfsGyk7aZmZWkqqqKli1bUlFRweZvRrWaRAQrV66kqqqKffbZp9btfHvczMxKsnbtWtq2beuEXQJJtG3bts53KZy0zcysZE7YpSvl2jlpm5lZo3f22Wczbty4T9XHl770JVatWlXndmPGjGH48OGf6tgbeU7bzMzqRakPr6pOxI/rtb9SRQQRwYQJE8odikfaZmaWTx988AHHH388hxxyCD169ODBBx+koqKCK664gn79+tGvXz+WLl26qf4zzzzDEUccwb777rvZqPvGG2/k0EMP5eCDD+bHP87+UKisrOSAAw7gwgsvpHfv3rzxxhtUVFTwzjvvAHDPPfdw8MEHc8ghh3DmmWcC8Pjjj3PYYYfRq1cvBg8ezPLly+v9nJ20zcwsl5588kk6dOjAvHnzWLhwIccddxwArVq1Yvr06QwfPpyLL754U/1ly5bx3HPPMX78eK688koAJk2axJIlS5g+fTpz585l1qxZPPPMMwAsXryYs846izlz5tC1a9dN/SxatIiRI0cyefJk5s2bx8033wzAF77wBaZNm8acOXM47bTT+NnPflbv5+zb42ZmlksHHXQQl112GVdccQUnnHACRx11FACnn376pp+XXHLJpvonn3wyTZo0oXv37ptGwZMmTWLSpEn06tULgDVr1rBkyRK6dOlC165d6d+//1bHnTx5Mqeccgrt2rUDYK+99gKyr8ANHTqUZcuW8dFHH9Xpq1y15aRtZma51K1bN2bNmsWECRO46qqrOPbYY4HNV2UXbu+6666btiNi08+rrrqK8847b7O+Kysr2WOPPYoeNyKKrvz+7ne/y6WXXsqJJ57I1KlTGTFiRMnnVp0ab49Lai5puqR5khYprTSQdIik5yUtkPS4pFYFba6StFTSYklDCsr7pPpLJd0if1fAzMxK9NZbb7H77rtzxhlncNlllzF79mwAHnzwwU0/Dz/88G32MWTIEO666y7WrFkDwJtvvsnbb7+9zTbHHHMMY8eOZeXKlQC8++67ALz33nt07NgRgLvvvrv0E9uG2oy01wGDImKNpGbAc5ImArcCl0XEXySdA1wO/Juk7sBpwIFAB+DPkrpFxAZgFDAMmAZMAI4DJtb7WZmZWaO3YMECLr/8cpo0aUKzZs0YNWoUp5xyCuvWreOwww7jk08+4f77799mH8ceeywvvvjipuTeokUL7r33Xpo2bVptmwMPPJCrr76ao48+mqZNm9KrVy/GjBnDiBEjOPXUU+nYsSP9+/fntddeq9fzBdDGWwS1qiztDjwHXAA8DbSOiJDUGXgqIrpLugogIn6S2jwFjAAqgSkR8flUfjowMCLO2+pABfr27RsN4X3aAytPLqndjn+fdu3/eZqZfRovvvgiBxxwQLnD2ExFRQUzZ87cNN/c0BW7hpJmRUTfYvVrtXpcUlNJc4G3gacj4gVgIXBiqnIqsDE7dQTeKGhelco6pu0ty4sdb5ikmZJmrlixojYhmpmZNXq1StoRsSEiegKdgH6SegDnABdJmgW0BD5K1YvNU8c2yosd746I6BsRfdu3b1+bEM3MzKisrMzNKLsUdfqedkSsAqYCx0XESxFxbET0Ae4HXknVqvjHqBuyRP9WKu9UpNzMzMxqoTarx9tLapO2dwMGAy9J+qdU1gT4EXB7avIn4DRJu0raB9gPmB4Ry4DVkvqnVeNnAY/V9wmZmZk1VrUZae8NTJE0H5hBNqc9Hjhd0svAS2Qj5t8BRMQiYCzwN+BJ4KK0chyyBWx3AkvJRuZeOW5mZlZLNX7lKyLmA72KlN8M3FxNm5HAyCLlM4EedQ/TzMzM/OxxMzPLpVWrVnHbbbcB2YNWTjnllDq1r4/Xde5ofoypmZnViy++Pq9e+5vS9ZBt7t+YtC+88EI6dOiQuwRcCo+0zcwsl6688kpeeeUVevbsyamnnkqPHtns64YNG7j88ss3vW7zN7/5DZA9M3z48OF0796d448/vsbHlTZEHmmbmVku3XDDDSxcuJC5c+dSWVnJCSecAMDo0aNp3bo1M2bMYN26dRx55JEce+yxzJkzh8WLF7NgwQKWL19O9+7dOeecc8p8FnXjpG1mZo3KpEmTmD9//qbb5e+99x5LlizhmWee4fTTT6dp06Z06NCBQYMGlTnSunPSNjOzRiUiuPXWWxkyZMhm5RMmTCj6Ss088Zy2mZnlUsuWLVm9evVW5UOGDGHUqFF8/PHHALz88st88MEHDBgwgAceeIANGzawbNkypkyZsqND/tQ80jYzs1xq27YtRx55JD169NjsTVnf+c53qKyspHfv3kQE7du359FHH+UrX/kKkydP5qCDDqJbt24cffTRZYy+NE7aZmZWL2r6itb2cN99921V1qRJE66//nquv/76rfb96le/2hFhbTe+PW5mZpYTTtpmZmY54aRtZmaWE07aZmZmOeGkbWZmlhNO2mZmZjnhr3yZmVmjMGLECFq0aMH777/PgAEDGDx4cLlDqndO2mZmVi/Ove3deu3vtxfuVVK7a6+9tl7jaEh8e9zMzHJr5MiR7L///gwePJjFixcDcPbZZ296Wci1117LoYceSo8ePRg2bBgRAcCMGTM4+OCDOfzww7n88ss3vdazoXPSNjOzXJo1axYPPPAAc+bM4eGHH2bGjBlb1Rk+fDgzZsxg4cKFfPjhh4wfPx6Ab3/729x+++08//zzNG3adEeHXrIak7ak5pKmS5onaZGka1J5T0nTJM2VNFNSv1ReIenDVD5X0u0FffWRtEDSUkm3KO+vWzEzs7J59tln+cpXvsLuu+9Oq1atOPHEE7eqM2XKFA477DAOOuggJk+ezKJFi1i1ahWrV6/miCOOAOAb3/jGjg69ZLWZ014HDIqINZKaAc9JmghcC1wTERMlfQn4GTAwtXklInoW6WsUMAyYBkwAjgMmfrpTMDOzndW2xn5r167lwgsvZObMmXTu3JkRI0awdu3aTbfI86jGkXZk1qRfm6VPpE+rVN4aeGtb/UjaG2gVEc9HdsXuAU4uMW4zM9vJDRgwgEceeYQPP/yQ1atX8/jjj2+2f+3atQC0a9eONWvWbJrn3nPPPWnZsiXTpk0D4IEHHtixgX8KtVo9LqkpMAv4HPDriHhB0sXAU5JuIkv+RxQ02UfSHOB94EcR8SzQEagqqFOVyoodbxjZiJwuXbrU6YTMzGzn0Lt3b4YOHUrPnj3p2rUrRx111Gb727Rpw7nnnstBBx1ERUUFhx566KZ9o0eP5txzz2WPPfZg4MCBtG7dekeHX5JaJe2I2AD0lNQGeERSD7KkeklEPCTp68BoYDCwDOgSESsl9QEelXQgUOweRtF7FBFxB3AHQN++ffN7H8PMbCdS6le0Po2rr76aq6++utr91113Hdddd91W5QceeCDz588H4IYbbqBv377bLcb6VKfvaUfEKklTyeaivwV8P+36I3BnqrOObB6ciJgl6RWgG9nIulNBd52o4Za6mZnZ9vDEE0/wk5/8hPXr19O1a1fGjBlT7pBqpcakLak98HFK2LuRjaZ/SpZwjwamAoOAJQX1342IDZL2BfYDXo2IdyWtltQfeAE4C7h1O5yTmZnZNg0dOpShQ4eWO4w6q81Ie2/g7jSv3QQYGxHjJa0Cbpa0C7CWNAcNDACulbQe2ACcHxEbH5NzATAG2I1s1bhXjpuZmdVSjUk7IuYDvYqUPwf0KVL+EPBQNX3NBPLx2BkzM7MGxk9EMzMzywknbTMzs5xw0jYzs9y65ZZbOOCAA/jmN79ZbZ0WLVoAUFlZmZsXg1THr+Y0M7P68fN6fp3ED2p+TMdtt93GxIkT2Weffer32A2UR9pmZpZL559/Pq+++ionnngirVu35qabbtq0r0ePHlRWVpYvuO3ESdvMzHLp9ttvp0OHDkyZMoVLLrmk3OHsEE7aZmZmOeGkbWZmubfLLrvwySefbPp94xu+GhsnbTMzy72Kigpmz54NwOzZs3nttdfKHNH24aRtZma597WvfY13332Xnj17MmrUKLp161bukLYLf+XLzMzqRy2+olXfCleIT5o0qWidNWvWANlofOHChTsirO3GI20zM7OccNI2MzPLCSdtMzOznHDSNjOzkkXs+HnsxqKUa+ekbWZmJWnevDkrV6504i5BRLBy5UqaN29ep3ZePW5mZiXp1KkTVVVVrFixotyh5FLz5s3p1KlTndo4aZuZWUmaNWu207xdq6Go8fa4pOaSpkuaJ2mRpGtSeU9J0yTNlTRTUr+CNldJWippsaQhBeV9JC1I+26RVM/vcTMzM2u8ajOnvQ4YFBGHAD2B4yT1B34GXBMRPYF/T78jqTtwGnAgcBxwm6Smqa9RwDBgv/Q5rt7OxMzMrJGrMWlHZk36tVn6RPq0SuWtgbfS9knAAxGxLiJeA5YC/STtDbSKiOcjW7VwD3ByvZ2JmZlZI1erOe00Up4FfA74dUS8IOli4ClJN5El/yNS9Y7AtILmVans47S9ZbmZmZnVQq2+8hURG9Jt8E5ko+YewAXAJRHRGbgEGJ2qF5unjm2Ub0XSsDRPPtOrEs3MzDJ1+p52RKwCppLNRX8LeDjt+iOwcSFaFdC5oFknslvnVWl7y/Jix7kjIvpGRN/27dvXJUQzM7NGqzarx9tLapO2dwMGAy+RJdyjU7VBwJK0/SfgNEm7StqHbMHZ9IhYBqyW1D+tGj8LeKw+T8bMzKwxq82c9t7A3WleuwkwNiLGS1oF3CxpF2At2apwImKRpLHA34D1wEURsSH1dQEwBtgNmJg+ZmZmVgs1Ju2ImA/0KlL+HNCnmjYjgZFFymcCPeoeppmZmfnZ42ZmZjnhpG1mZpYTTtpmZmY54aRtZmaWE07aZmZmOeGkbWZmlhNO2mZmZjnhpG1mZpYTTtpmZmY54aRtZmaWE07aZmZmOeGkbWZmlhNO2mZmZjnhpG1mZpYTTtpmZmY54aRtZmaWE07aZmZmOeGkbWZmlhNO2mZmZjlRY9KW1FzSdEnzJC2SdE0qf1DS3PSplDQ3lVdI+rBg3+0FffWRtEDSUkm3SNJ2OzMzM7NGZpda1FkHDIqINZKaAc9JmhgRQzdWkPRz4L2CNq9ERM8ifY0ChgHTgAnAccDEUoM3MzPbmdQ40o7MmvRrs/SJjfvTaPnrwP3b6kfS3kCriHg+IgK4Bzi5xLjNzMx2OrUZaSOpKTAL+Bzw64h4oWD3UcDyiFhSULaPpDnA+8CPIuJZoCNQVVCnKpUVO94wshE5Xbp0qeWpmJlZQ5RmVess4scltTv3tndLavfbC/cqqd2OVKukHREbgJ6S2gCPSOoREQvT7tPZfJS9DOgSESsl9QEelXQgUGz+OoqUERF3AHcA9O3bt2gdMzOzevXzEpdZ/WDHpalaJe2NImKVpKlkc9ELJe0CfBXoU1BnHdk8OBExS9IrQDeykXWngu46AW99qujNzMx2IrVZPd4+jbCRtBswGHgp7R4MvBQRVVvUb5q29wX2A16NiGXAakn90zz4WcBj9XkyZmZmjVltRtp7A3enRNwEGBsR49O+09h6AdoA4FpJ64ENwPkRsXGC4QJgDLAb2apxrxw3MzOrpRqTdkTMB3pVs+/sImUPAQ9VU38m0KNuIZqZmRn4iWhmZma54aRtZmaWE07aZmZmOeGkbWZmlhN1+p62mZlVLzdP/vqwbUntduRDRKw4j7TNzMxywknbzMwsJ5y0zczMcsJJ28zMLCectM3MzHLCSdvMzCwn/JUvMzNrkL74+ryS2n2OzvUcScPhkbaZmVlOOGmbmZnlhG+Pm1md+clfZuXhkbaZmVlOOGmbmZnlhJO2mZlZTnhO28yszPzVJqutGkfakppLmi5pnqRFSitQJD0oaW76VEqaW9DmKklLJS2WNKSgvI+kBWnfLZK0Xc7KzMysEarNSHsdMCgi1khqBjwnaWJEDN1YQdLPgffSdnfgNOBAoAPwZ0ndImIDMAoYBkwDJgDHARPr84TMzMwaqxpH2pFZk35tlj6bvg+RRstfB+5PRScBD0TEuoh4DVgK9JO0N9AqIp6PiADuAU6utzMxMzNr5Gq1EE1S03T7+23g6Yh4oWD3UcDyiFiSfu8IvFGwvyqVdUzbW5YXO94wSTMlzVyxYkWtTsTMzKyxq1XSjogNEdET6EQ2au5RsPt0/jHKBig2Tx3bKC92vDsiom9E9G3fvn1tQjQzM2v06rR6PCJWSZpKNhe9UNIuwFeBPgXVqmCzJY2dgLdSeaci5Wa2k/AqabNPpzarx9tLapO2dwMGAy+l3YOBlyKi8Lb3n4DTJO0qaR9gP2B6RCwDVkvqn+bBzwIeq79TMTMza9xqc3t8b2CKpPnADLI57fFp32lsfmuciFgEjAX+BjwJXJRWjgNcANxJtjjtFbxy3HaQgQMHljsEM7NPrcbb4xExH+hVzb6zqykfCYwsUj4T6LF1C7P69+yzz7J+/Xq++MUvbiqbPHkyzZo146ijjqpVH34xhpk1JH4imjVaXbt25dJLL2XcuHGsXr2aiy66iOXLl/OLX/yi3KGZmZXEzx63RqtLly6MGzeOVq1aMXv2bFq1asW4cePo0qVLuUMzMyuJk7Y1WlVVVQwdOpRVq1bRu3dvVq1axdChQ6mqqqq5sZlZA+SkbY3Wq6++ynnnnceoUaNo2bIlo0aN4rzzzuO1114rd2hmZiXxnLY1WgMGDNiqbNCgQTvk2P4+spltD07a9qmUmpymdD2kniPZtqlTp+7Q45mZbQ++PW5mZpYTTtpmZmY54aRtZmaWE07aZmZmOeGkbWZmlhNePW5A6c/YHlh5cknt/IxtM7O680jbzMwsJxTRsEcgklYAr5c7ju2kHfBOuYNogHxdivN1Kc7XpThfl+LycF26RkT7YjsafNJuzCTNjIi+5Y6jofF1Kc7XpThfl+J8XYrL+3Xx7XEzM7OccNI2MzPLCSft8rqj3AE0UL4uxfm6FOfrUpyvS3G5vi6e0zYzM8sJj7TNzMxywkm7DCTdJeltSQvLHUtDIqmzpCmSXpS0SNL3yx1TuUlqLmm6pHnpmpT2FJxGSlJTSXMkjS93LA2FpEpJCyTNlTSz3PE0FJLaSBon6aX0/5jDyx1TKXx7vAwkDQDWAPdERI9yx9NQSNob2DsiZktqCcwCTo6Iv5U5tLKRJGCPiFgjqRnwHPD9iJhW5tAaBEmXAn2BVhFxQrnjaQgkVQJ9I6Khfxd5h5J0N/BsRNwp6TPA7hGxqsxh1ZlH2mUQEc8ApT3HsxGLiGURMTttrwZeBDqWN6ryisya9Guz9PFf2oCkTsDxwJ3ljsUaNkmtgAHAaICI+CiPCRuctK2BklQB9AJeKHMoZZduAc8F3gaejoid/pokvwT+FfikzHE0NAFMkjRL0rByB9NA7AusAH6XplPulLRHuYMqhZO2NTiSWgAPARdHxPvljqfcImJDRPQEOgH9JO30UyqSTgDejohZ5Y6lAToyInoD/wJclKbjdna7AL2BURHRC/gAuLK8IZXGSdsalDRv+xDwh4h4uNzxNCTpdt5U4LjyRtIgHAmcmOZvHwAGSbq3vCE1DBHxVvr5NvAI0K+8ETUIVUBVwV2qcWRJPHectK3BSIuuRgMvRsQvyh1PQyCpvaQ2aXs3YDDwUlmDagAi4qqI6BQRFcBpwOSIOKPMYZWdpD3SIk7S7d9jgZ3+WyoR8d/AG5L2T0XHALlc4Or3aZeBpPuBgUA7SVXAjyNidHmjahCOBM4EFqQ5XIAfRsSE8oVUdnsDd0tqSvZH9tiI8NebrDqfBR7J/v5lF+C+iHiyvCE1GN8F/pBWjr8KfLvM8ZTEX/kyMzPLCd8eNzMzywknbTMzs5xw0jYzM8sJJ20zM7OccNI2MzPLCSdtszqSFJJ+XvD7ZZJG1FPfYySdUh991XCcU9ObjqZs72OVm6QfljsGs/ripG1Wd+uAr0pqV+5ACqXvctfW/wEujIgvbq94GhAnbWs0nLTN6m49cAdwyZY7thwpS1qTfg6U9BdJYyW9LOkGSd9M78peIOmfC7oZLOnZVO+E1L6ppBslzZA0X9J5Bf1OkXQfsKBIPKen/hdK+mkq+3fgC8Dtkm4s0uZfU5t5km5IZT0lTUvHfkTSnql8qqT/K+mZNHI/VNLDkpZIui7VqUjvML47tR8nafe075j0AocFyt4zv2sqr5R0jaTZad/nU/keqd6M1O6kVH52Ou6T6dg/S+U3ALspe7f0H1L7J9K5LZQ0tA7/3M3KLyL88cefOnzI3oXeCqgEWgOXASPSvjHAKYV108+BwCqyJ5ztCrwJXJP2fR/4ZUH7J8n+oN6P7JnJzYFhwI9SnV2BmcA+qd8PgH2KxNkB+H9Ae7KnY00mez85ZM8w71ukzb8AfyV71zDAXunnfODotH1tQbxTgZ8WnMdbBedYBbQFKsjePHVkqndXumbNgTeAbqn8HrKXxJCu7XfT9oXAnWn7euCMtN0GeBnYAzib7ClXrVO/rwOdC/8ZpO2vAb8t+L11uf998sefunw80jYrQWRvH7sH+F4dms2I7J3h64BXgEmpfAFZYttobER8EhFLyBLR58meIX1WerzrC2TJcL9Uf3pEvFbkeIcCUyNiRUSsB/5A9k7hbRkM/C4i/p7O811JrYE2EfGXVOfuLfr5U8F5LCo4x1eBzmnfGxHxX2n7XrKR/v7AaxHxcjX9bnxhzCz+cX2OBa5M12EqWYLukvb9Z0S8FxFryZ4r3bXI+S0gu5PxU0lHRcR7NVwPswbFzx43K90vgdnA7wrK1pOmndILUD5TsG9dwfYnBb9/wub/LW75bOEARDbyfKpwh6SBZCPtYlRD/NW1qeuzjQvPY8tz3Hhe1Z1TbfrdUNCPgK9FxOLCipIO2+LYhW3+cdCIlyX1Ab4E/ETSpIi4toY4zBoMj7TNShQR7wJjyRZ1bVQJ9EnbJwHNSuj6VElN0jz3vsBi4CngAmWvLkVSt/QWp215AThaUru0SO104C81tJkEnFMw57xXGo3+j6SjUp0za9HPlrpIOjxtnw48R/a2sgpJn6tDv08B301/ECGpVy2O/XHBdesA/D0i7gVuIqevZ7Sdl0faZp/Oz4HhBb//FnhM0nTgP6l+FLwti8mS12eB8yNiraQ7yW4Rz04JawVw8rY6iYhlkq4CppCNUCdExGM1tHlSUk9gpqSPgAlkq6+/RbZwbXdKe0PSi8C3JP0GWAKMSuf1beCPknYBZgC319DPf5Dd4ZifrkMlcEINbe5I9WeTTWncKOkT4GPggjqeh1lZ+S1fZrZdSaoAxkdEj3LHYpZ3vj1uZmaWEx5pm5mZ5YRH2mZmZjnhpG1mZpYTTtpmZmY54aRtZmaWE07aZmZmOeGkbWZmlhP/H8Cbjh3RfvVqAAAAAElFTkSuQmCC",
            "text/plain": [
              "<Figure size 576x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "bic = np.array(bic)\n",
        "color_iter = itertools.cycle([\"navy\", \"turquoise\", \"cornflowerblue\", \"darkorange\"])\n",
        "clf = best_gmm\n",
        "bars = []\n",
        "\n",
        "# Plot the BIC scores\n",
        "plt.figure(figsize=(8, 6))\n",
        "spl = plt.subplot(2, 1, 1)\n",
        "for i, (cv_type, color) in enumerate(zip(cv_types, color_iter)):\n",
        "    xpos = np.array(n_components_range) + 0.2 * (i - 2)\n",
        "    bars.append(\n",
        "        plt.bar(\n",
        "            xpos,\n",
        "            bic[i * len(n_components_range) : (i + 1) * len(n_components_range)],\n",
        "            width=0.2,\n",
        "            color=color,\n",
        "        )\n",
        "    )\n",
        "plt.xticks(n_components_range)\n",
        "plt.ylim([bic.min() * 1.01 - 0.01 * bic.max(), bic.max()])\n",
        "plt.title(\"BIC score per model\")\n",
        "xpos = (\n",
        "    np.mod(bic.argmin(), len(n_components_range))\n",
        "    + 0.65\n",
        "    + 0.2 * np.floor(bic.argmin() / len(n_components_range))\n",
        ")\n",
        "plt.text(xpos, bic.min() * 0.97 + 0.03 * bic.max(), \"*\", fontsize=14)\n",
        "spl.set_xlabel(\"Number of components\")\n",
        "spl.legend([b[0] for b in bars], cv_types)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAACLCAYAAABflUO5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAtyUlEQVR4nO2deZxU1Zn3v6f2pfduloamGxAUEIXQIh1eMUaNIjFxjA5jxjcTmZkXQ/SNxmQmM68z72tmMomZNxqcMb7BiYOJISFMiHE3bsGNEBFUICxBloYGgd6rurZbt+55/7gLVdXV3dUK9ML5fj79qap7zj333Nv3/uqp5zznOUJKiUKhUCjOPK6h7oBCoVCcrSgBVigUiiFCCbBCoVAMEUqAFQqFYohQAqxQKBRDhBJghUKhGCKUAPeDEEIKIaYNcR9uEUK8MZR9GCzCZLUQolMI8VYR9R8VQnzLen+ZEKLl9PeyYD82CCH+usi6Q35vKEY+o16AhRCXCCE2CiG6hRAdQog3hRDzz+Dxi36oP2T7PiHE/xZC7BFCxIQQR4QQzwkhrsqqc1AIoQkhavL2fdcSksnW50etz5/Nq7fS2n5Lkd26BPgUUCelvPgjneAIRwgxVgjxcyHEUesefFMIsWCo+zUcEULcI4T46VD340wyqgVYCFEGPA38O1AFTAS+CaSGsl+nmF8C1wF/AVQCU4AHgE/n1TsAfN7+IIS4AAgWaO+PwBez6nmAPwX2DaJPDcBBKWVsEPuMVkqAzUAj5j34Y+AZIUTJkPZKMTyQUo7aP+AioGuAOn8J7AI6gd8ADVllEphmvfcD3wMOAceBHwLBrLrXAe8CEUyxWgz8C5ABkkAP8KBVdwbwItAB7AGWZrVTDTxptfMW8M/AG330/UoggWlp9neOB4F/ADZnbfsecLd1jpOtbY9a248Blda2a4HngDeAW4q45n9lnW/GOudvArfkn0PetX0U+Jb1/jKgpZ/2JfBlYC8Qta7POcDvrGu2DvBl1f8fwPvWtX4SmJBV9ilgN9ANPAi8Cvz1YO+ND3FfRoDGPsrcwP+y7qEosAWYZJUtxBTzbut1YdZ+G4BvARut6/6UdS+tsY632f4/Z/X/K8B+oA34v4DLKnNZ90szcAL4CVBulU229v0i5rPQBtyd1a4L+Dur/+3W/6NqoH0xnxcNSFv9f8/afovVxyimEXHzUOvKqfwb8g6c1pODMusm+DFwDZaoZJX/ifVwzgQ81k23Me8mtUVipfUAVwGl1g3+HavsYuuh+JR1A04EZmQ9GNkPdRg4DCyzjjnPuhHPt8rXWjdtGJgNHKFvAb4X2FDEdTiIKdZ7rHN1W31ooLcAfwt4GFhhbVuHaTk7AgzUA11AfR/HuyW7z/mfC1zbRxmcAD9p/W/Px/w18zIwFSgHdgJftOpebl3beZhfoP8OvGaV1WAK042AF/gqoNv/q8HcG4O8J+difkGV91H+N8B24DxAAHMwhbQK84vgC1Z/Pm99rs66z97H/DKyr8Mfrf+7B1NEV+f1/7dWu/VWXfvc/9JqayqmBf8r4DGrbLK1739g/oKaY/0PZlrldwKbgDrrmq8Cfl7kvvcAP817ViLAedbnWqznZLT8DXkHTvsJmg/Qo0CL9YA9CYyzyp4D/iqrrguIY1k69kNmPQgx4Jysuh8HDljvVwHf7+P4G8gV4D8DXs+rswr4P5jCmMYSb6vs2/QtwD8C1mZ9rsIUxm4gmbX9oPUg/gPwHUxr40XrwSwkwJdgWpTlmNZ+kCItYKudWzi9Avzfsj5vAb6R9fk+YKX1/hHgX7PKSqzrOxnTZbMpq0xY94gtQkXdG4O8F8swxfXv+6mzB7iuwPYvAG/lbfsdJ78UN5Brid4HPJf1+TPAu3nXcXHW5y8DL1vvXwa+nFV2nnXdPJwU0bqs8reAm6z3u4ArsspqB7HvPfQW4C7gBrJ+bY6mv1HtAwaQUu6SUt4ipazDtCgnYFqzYFqADwghuoQQXZg/UwWmBZvNGCAEbMmq+7y1HWASxftIG4AFdjtWWzcD4632PJjWqU1zP221Y97g9rl2SCkrMP2N/gL1HwP+HFMQf9JXo1LKN6y+/APwtJQyMeBZnVmOZ71PFPhs+1cnkHX9pJQ9mNdsolV2OKtMknvdi703ikIIEcT81bRJSvmdfqr2dS/lnItFc15/ir0uNvn32YQ+jtWMeV+Oy9p2LOt9PKvtBuDxrOu2C9MdVcy+OUhzDOHPgC8BHwghnhFCzChUd6Qy6gU4Gynlbkxra7a16TBwq5SyIusvKKXcmLdrG+YNfH5WvXIpZUlWO+f0ddi8z4eBV/OOWSKlXAG0Ylrpk7Lq1/dzSi8D84UQdf3UOdkRKZsx/WhLMH9W9sdPga/Rj1APghjmFxgAQojxp6DNYjiKKQj2ccOYP+ePAB+QdZ2FEILc617svTEgQgg/8GvruLcOUL2veynnXCzqrTY/LPn32dE+jlWPeV9mC3pfHAauybtuASllMf3Mf1aQUv5GSvkpTENjN6b7YtQwqgVYCDFDCPE1W6CEEJMwfWebrCo/BP5eCHG+VV4uhPjT/HaklAbmP/77QoixVt2JQoirrSqPAMuEEFcIIVxWmf1NfRzTl2bzNHCuEOILQgiv9TdfCDFTSpnBFMZ7hBAhIcQssiISCvTrBUw/3q+FEAuskDQv0NTPZfkr4HI5cITCv2H6tF8boF4xvAecL4SYK4QIYP7UPBP8DPP/MtcSwW8Dv5dSHgSesfr0OSvS4yuYv0Jsiro3rLJbhBAH+yjzYkaqJIC/sO6l/vgR8M9CiOlWPPWFQohq4FnM++bPhRAeIcSfAbMw76cPy98IISqt5+IO4BfW9p8DXxVCTLGiNb4N/EJKqRfR5g+BfxFCNAAIIcYIIa4rsj/HgclCCJe17zghxGetL84U5uBcpuizGwGMagHGHDldAPxeCBHDFN4dmJYdUsrHge8Ca4UQEavsmj7a+gbmwMQmq+5LmL4xpJRvYQ6qfR/T//oqJy2IB4AbhTkp4d+klFHgKuAmTEvjmNUH22VwO+ZPsmOY1vrqAc7xc5gP4U8x/WUHMF0aiwtVllLuk1K+PUCbtjvjZeuneQ5CiHohRI8Qoj/rPLutPwL/hHnN9mL6k087UsqXgX8E1mNavOdgXneklG2Y4XX3YrolpgNvZu07mHtjUva+eSzEjCS5CuiyrluPEGJRH/Xvxxz4fAFzAOoRTP9nu9XO16z+/i1wrXUeH5YnMH3o72J+IT1ibf9PTHfVa5j3UxL4n0W2+QDmOMsLQogo5jNXbNzzf1mv7UKIrZj69DXM56QD+ASmr3rUIAo8XwqFYhAIIV4A7pBS7hrqvhSLEEIC06WU7w91X85mPEPdAYVipCOlvGrgWgpFb0a7C0KhUCiGLcoFoVAoFEOEsoAVCoViiFACrFAoFEPEoAbhampq5OTJk09TVxQKhWJ0smXLljYp5Zj87YMS4MmTJ/P22wOGkCoUCoUiCyFEwZQCygWhUCgUQ4QSYIVCoRgizioBjhkZnoq2EzNG1XRyhUIxQjmrBPiVWBf3d7TwSqyr33pKqBUKxZngrJqKfHm4Iue1L2yhBvhMafVp7pVCoThbOasEOOxyFyWoxQq1QqFQfBRGvAvidLgLbKEOu9yDPlZfdZRbQ6FQ5DPiBbhYv+6p4LmeDu7vaOF7bYf7FNK++nMm+6lQKEYGI94FUay7IGZkeCXWxeXhil6WbbF1hfW6IdHNvFhpQXdGX/1Rbg2FQpHPiBfgfL9uX+I5mIG153s6eLDzKJo0uKHs5OzBxSVVgLlwVVOwlKei7b2O05efuVj/s0KhOHsY8QKcTczIsLK9hZfiXUCu0A7GApV5rzZhl5vFJVW8EuvitXg3D3YedY4zGAtboVAoYAT7gAsNar0S6+KleBdXhip6CW32wNpAA2LXlFRxV1Ud15RU9aprW9I7kjFur5zgHGcgH++ZHIRTA34KxchgxFrAhVwK2VZuoQgG20K1992W7OHO6rpedbPdBU9F23OOc3m4gm3JHl6Kd+ERwnFLXB6uICUNNGkQMzK92jxTscX9/QpQKBTDixErwIVcCv35WbMFMFtELwyU9CtS+ccJu9zcWV0H4Ijc8spaNiWiCODBzqP4hKtXm2dqEK6/XwEKhWJ4MWIFuNjBN5t863h5ZS1gDqb1RaE27W03lFazLRXjpXgXOpIN8W5uq5zAXVV1BYXvTA3C9fcrQKFQDC+GvQ94IH+mXW7H6L4S6yq4T9jlpilYysr2Flp1jU2JKC/Fu9iUiPZ5rEJ+Xfs493cc4UQmDcB0X9DxGReawFHMeRRbZyD6mkSiUCiGH8PaAh7In5ldfnuW9dmXv/WhzqNsiHejI/l69SQg1yVg75eWkoWhMhq8fm4sqcGP4P72Fmb6guxOxQGo9/hJSINql5f5gVLO8QVxCUFfFOMDVjkoFIqzi2EtwAP5M7PLF5dUOREOKWnkRCjYTPcF2RDvpsHrz3EtpKXBkbQGSBr9JexMxXhfM4XWAJ7oaWOnlmCz20uFy81MX5BDeooWXaMFjQc6jjDNF6Te62eKN0Ct109YuNgQ73aOUYwPWE3WUCjOLoZEgIuNmR3In1mo/JVYFz/oPMpdVb2jG64rraHU5UGTBvd3tNCeSVPt9rJHi5OREgPzrzuTZro/hE+YHpoqt4eEYXBAT3Esk6YpUIomDQAmeXzM9IfYryXYpyXY5fHiES5O6Gm2pnpIGBk+XVpd1PmqyRoKxdnFkAhwoZ/ahUR5IEEqVN6fFRl2ubm2pIpdWpyPB8o4nE7RkUlT4fJgINkUj3BATwEghMCNoM7rpyWdojFYSqXmoT2jM87j5d1kzOn3fi3B1pT5ucRVyhSfj2O6xjRvgMejbWxL9vBmMppzvgqFQjEkAlxIJE+V/7O/6Ij2TJofdnyAVwjiRoYGrw8hBPu0BBkkB/QUkzw+XAgyUrI51cMUXeOAnqIJ8AkXzXoKVxKa9RRh4aLDyCDTSeo8PgRQ5/VzQEvybipGtctDu6EDUOv2gZS062k2JiK9rGE1k06hOPs4IwKcLy6DtVw/CrawH0mn+EMqxg4tzmS3j4MZjYihU+PxsjnZw/xACfP8YQ6mU3QYOmM9XpoCpYzzeCFpCitAg65R5nYzzx0mg+RQOkWnkaHTilwYn06SkZIGj5/Z/hA7U3Eq3R62pmI8H+vglXgX71rWcrb1ryZPKBRnH2dEgIuxbj+M/zNmZHi+pwOJOX24kOV4aaic7ckYJzIa5/qDlLjc1Hn9iESUA3qKGkto67x+fhvrosPQCQnTAp7mD7InFeeAnqJS8zgWcLMODR4/zXqKBo+fTiNDvcfHOI8PwHFH1Hp8XBquQJMGPuFiii9AwsigSckYtwcpJUIINXlCoThLOSMCfLqs2+d6OviBlRDHX2D2WUZKXo13IYSkzu3HJQTn+c1TbgqVMU5LMsUXwCdc7EnFaTd0wsJFTBpsTcXwCRftGdOF0J7RWRgq4wNdo1lPUeF248LPhYEwtbrPaSdmZDihp6l2e5jiCwCm6+I8fwhNGhxKp4hmdJ7t6SApJZeGyouaQt0ULGVTIqpcFArFKOKMCPBHHd0fKD/vZVkils0b8W52pGLUuk1fL4AmDQ5oScelYGOLZZ3Xz34tQXtGp87rp87rx5WAxmApcSND1Mgw2xekK5OhWU9hJKHa7eF9LcE0X5CWdIpmPUWtx+dEUWjS4H0twbG0xuGMZvZdgzdd3XgQdBp6wdl22THNV4YqerkolN9YoRjZDOs4YJu+oiYkcHvlBC4NlTtCBGY+3/ZMmmO6xkSP3xFfgPe1BJuTPTRYlizAeVbImW2ltmd0mvUUY9NeZvnDXGq1+1TUdFH0aDoa4ENYLgmznf1akktCZY5LY6flishI6bglbCo9Hsa6vayJHGe3lsjJPZyftwLMKdMXBkpOy8ClQqEYGkaEAOe7MLIHre6qqmNTIuoIEeDk6b3AFyJhGI57QJMGJ3Rz+nC120Otx3Qd2FbxFF+AA1rSEdR8Jnl9dKRM8bVdFQCVlvXZbui8EY9wVUklB7Qkm5M9ZrnobZ26EXiFC69lx6elpFXXeLjzA/57+VhnVp89iaOQpdsULOXKUEW/+SwUCsXwZUQIcL4Lo69Bq6ZgKa/Fu7mhtJodqTh+Idhkxd+e5w854lrt8jDVF3TEbGcqxuZkD5o0cAvBPH8YtxBM8wVzXBZuhBluJuHCQJgdqTiGlBzOaMzzhxHpFO2GzsZ4hPnBUuYjOaGncwQ9iCCBpDWTZmcqxnn+EBlgrNvL/249yG4tAcCd1XW93BCQa+na+SwGyuimUCiGJyNCgPMpNGj1mdJqnoq282DnUS4NljHTF8IjhOPLBdPPe9yK692vJejM6DRmWY+tmTSHdY35gRJm+cMA7EnF2ZSMMklPcVjXnLpjdC+1Hh/jPF48yRi1Hh8H06bQZrsvJnjStMbSxKVBg8fP/GApW6wIjMO6hhvBBf4wz/R0sFtLMMMXZHllbU7O4uWVtQWzrKmpywrFyGZECrBtEdvZw7J/qksp2a3FCbvc7NMSHNBTjEv7HD9vU6iMGi3BPi1Jh6GjxSUJmWG2L0TEmjSRjT04d8SyYr1AmpNibbsiuo0MHYbuWLgZaS5otC0ZIy4NwsLF3ECY5nSSSreHSrcHtxCOeyTkcjHV6+fvqicxxuMrKmfxqZi6rAbyFIqhY9ino+yP/HSRYZebT4QrkIDHEremQClTfAFiRobXYl2kLWG04317DJ0OazrxIV2jyuUhI6WT6yEtJcd1jZm+EGHhIm0dWzckbiAmDapdHiZ4vABUus3vtAyS12JdzPKHqHZ5iEmDbUnT1WGHuM3yh51IiX1agv3pFN9uP8wvI62A6YboK79wNh8ljeVASykpFIrTx4i0gG0KLQPUnknTkk5R4nLTkk5R5nLzZLSdcpebo5k0JKDGEsu4lExwe+jOGEz1BYkaZmhZR0pnj5bgynAF25IxDugpIkaGmDQICYEXwQeGKcU+oM7ro103recxHi/13oATLww4g3J1Xj81afPYtmUNZphaqduDC/ijluCPWsKJaz7VKz4XuobZrwqF4swxogUYTB/tS/EuZxmgtxNRdmhxYkaGA3oKHwINiZYxfbCNwVK8QjiDY+M9Puq9Aeq8fjZbydlty/aNeIRLQmVE4hkaAyW8nugmLiUgqbAiG0pcbt6zcgQDdGUyzAqFTb9zAsrcbt7XEjR4AzSnkwBM8wVzYoTtJEBXhSu43hdCkiuItvBq0ii4EnNTsLQoS7kQKgObQjF0jCgXRKEVK16KdzHDF3RCsa4qqWSWL8S8QAlNgVIuC5XjQ5DGnBpsW45jPV7mB0o4zx/iPH/ImUBR7fKQwQwzuyRUxlFdo90wLeKE5b7wAl0yQ7nbQ0tGI4RgljdIlctDs55iUzyCVwhqPT7eS8XZnOxhcyJqviZ72JOKsydlfknY4lvmctEUKOOGsjFcU1LlrOxhn+f9HS1I4LbKCY7Fb2/flIgWXPFZrY6sUAxvRpQFnP9TO3ug6tV4N37h4vJwBef7w7gwQ88AbiirceJ8sy3OKR6/E2qWQTLPinxowI9bmHG6dtywISWzfaZQB12CDzI6FW43UcNNh5HhWCZNh6FT7fJwQE9hxCNUuN2Md3vozmRIZomgPdHDzrRmZ007ZLks8idiaNLgtsoJjjDf39GCzzpXu06hawSoiRoKxTBmRAlwXysUXxgocZKsA3wsEOa1eDcBl69XGwe0pCN6B/QUNVqCNj3tCHK2INqRElUuD4czGi4hKHd7HN+uG0GJyxTg8R4P6DDW40HqOEl7bBJWhIUdilab9lHn9TMu7aPS5ebNZJQFgZJe5/l8TwcPdh7l9soJvVbWKDarnPLvKhTDkxElwIUEJzskDSAlDaZ4A7xONxkpcQvBAS3pTMjIzvnQkk45eYCrXR4uDITRk5LDukaFcNNhieZkr5/SjNsR3kqXmyles51DusYkj49juk6HodOhnRRae5KGF5juDRJwuSzLWjjW+Xn+EC/2dBIxMqyPtnN3oCTnPKV1nvbrYJPUK8tXoRi+jCgBLkTMyPBcTwcCU6R+0HmUPak4FwVKeSsZJZJJk8bMz1vn9TtpIe2ohOZ00kkteVTXiBlm+JmwGqy3UkxKRwKh08hwnnCRsba5EFZYm4vJXh8xQzLfGuxzfLxWeJr9ReDGDJPzIJjuC3IOAZZX1uac1yuxLj4RKndcKwqFYnQxrAW4mEkC9hpwYCbmsbOGTfcFOagl2J+d1yFhDr4BbE72UK0laDd05gdKqPX4yCDpsFJSzg+UsldLOMnU4eQkjArhJoOkwbKCM9KceNGS0ejOGBzJaERjOpeHK3PSXtpkkI4QV7g9fDxUxtUlVTnnnZ3rQlmxCsXoZFgLcDFrxzUFS7ksWM65/iCLS6pYXFLFhYESUtJgv56i1u2l2u2h24rxbdZTzA+UOP5eeyDOTtZj+4PfTESISQMDmOcPOzPfAEpdbifRjl2/ymVeyk4rPrjDyLAlEeXScEWOxQ3Q4A3gRlDh9lDu8tCZ0VkfaXVWdj4VCdrVDDeFYvgzrAW4v1H+bcke7qw2M6FtSHQzL1jqCI09eHVb5QTmBsI8He1gr7XMfL3HxzRfkGm+YI5luicVZ4ovQGOwlBM9aWLW9OFmPYULGOP2cljXqBRuqjwexntN14Qt4hcGwrwRj9Bu6NR7fAgEjUFz9eQ34t0c1jUnn0QT5oSNgHBR5fbw/7o+AHBimZuCZvhcndfHcz0dfa720R8qVaVCMfwZ1gLc1yh/do6EvkT6QWtp+mm+EDeWufm74/sBqHZ5HGs3g+R9K/uYnQ2tOZ1yphdfEipzZsLVeLyOr7gzFXemONu+XJ9w8clwBVsSURqzvgz2pOKO5SyBi/3mIFtIuKl0e1gQLGVXKs65/qBzDq/Gu9kQ73bOp9BqHwOhZrgpFMOfYS3AhcgOPSs2FOudZA9HrJUoIjJDq67x+0TUWbF4fqCE+YGSnLCzqb4AIZfb8eHWef1OTHCdx8eRdApNGk6SH4DmdDJn/bgpvgBTfAGOpFMczmhUujyE3W48QjiW7/taopcFb6ePXxQs65WEfTDXSVm+CsXwZsQJ8EAU8n02BUsZ6/ZwIqNT6fayLRmj3dApw0Wlx+tME+4wdCpdbsIuF5uTPbgxw8XsBTvbDZ0pHj8G0JwxRdVeSSObA+kknUaGDJJZ/jAXBUsptaYjl7nc/Czaym2VE7irqq7gSheLS6qciRbKf6tQjF7OuACfisGh/vybhco2JaKcyOjM8AX52+pJ/Kang0e6jyOEuaRQOm7gsqLMyqzpxA0ePxnMrGhbLGu52uWhKVRGWkqMhLlD9tpy03xBZ1AOIGkYtKRTtGfSfKF8LPOCpaSlZLzX3yuX8am8PgqFYmRwxgX4VAwOFfJv9pWYJmZkcqbyhl1uri8bQ3M6xUvxLiZ7/By0BLPW7WOqL0CF202nrp9cUsjtQUcyxm2GsIVdbsZ7fGxO9jA+nWSWP0zKyPC+lmSaL4jfZbokLgqWkjAMXo53cXmmEp9w4RP9n7caPFMozh7OuACfisGhbP+mLbxRQ+c/uo5xW+UEbrQWt4STS9ffVjkBgKei7TQFS5lhJeH5RKic1+LdJKXBhf4w7Zk0z/d0cigTp1y4iWQy7EknqHX72JqKkZKSeq+fiDXzLmpkOKFrHEqn+IMW5/JQBXdW1VHv9RO0kuNM8PppCpbmJI8/nddHoVCMDM64AJ/qwSHbYrwsVA6cHMCysT+npcH32g6zIdHNJcEy3khEuCxUzjUlVc5qxGAKentGxycEryciLAmWcGEgTHcmTUXazeJwJSVuNz4Es3whPhGuYJzHhxvT1ZEvsPb5PhVtHzCmObt+/mofCoVi9DHiB+Gyl22fFyjtZTnaA1opabAh0Z1TtiHe7eyTvQDmg5bFPD9Y5mz7r2ibeTwk0jB4ItbBXVV1zLES6EBvl0G2wGZbtvb2lDScWXz5+9pJeLKXq1coFKOLES/A2RZ1oeXb85P1CODSUDlz4iWkrdU0bLGz27Bfsyd2bE1G2RDvRjCwm6AvgbX7aVvDt1uRENmibB/XzjyxIxVjsTH4iRgKhWL4M6ISsg9Ef+ubhV1ubiwbww1lYxjj8XFj2RhKXB4e7DyKBCckzBZBgF9GWllvrc/29epJ3FVV50wXzhbT/ITndj+E1W6h1YzttuxE6vl9v6akiitDFWyId6v12hSKUcqwtoCLCcnq62d+MeRbu9l+WsCxXn0F1meD3onT7SiMlDS4vXKCI9b5/txiJo/kTzhRKBSjj2EtwMWEZOXXGcwAX74Q5otgSho5Lof8L4Ts+nY/7Gxsd1XVOeJrZzbr7zz6y3WsUChGJ8NagIuxaAdr9fZnVeeHt/nzZqPZIttj6OzXkiyvrM3xPwO9ZradisxmCoVidDKsBbgYC3CwVmKxEx0K1bMFdKuVgQ3g7jENvfrRl1VdyI2iZr4pFGcvw1qATwfFWsyF6tki2xQsxdMpclaw6IuBviDyhV4JskJx9nDWCXCxFnN/9cZ4fI7l+1HJF3o1FVmhOHs46wR4uDHQQKBCoRi9KAEeZqjIB4Xi7GFUTcQYaux435iRyXmvUCgUhVAW8Ckk238LKF+uQqHoFyXAp5BC/lvly1UoFH2hXBCnENt/mz3dWIWSKRSnlkgkxcMPbyESSY344ykBVigUI4q1a3dw661Ps3btjqL3iURSPPDA73nggU2DFtIPc7xiUS4IhUIxbIlEUqxe/S7JZJpAwMuyZXO56abZAM5roX3Wrt3hlK9du4NEIs2dd/4GgGDQy/LljUX3YaDjfRSUACsUimFJJJLiy19+hjVrtjvbgkEPy5c3OgKaLbZlZf6C+9x669OsXLmYlSsXA3LQQlpW5h+UYA8GJcAKheKMYlq17wCCZcvmUlbmz9r+LiBZtuxjrF27gzVrttPYWMtnPnMugYCHRCJNJJJy9rHdAwDLlzeyevW7rFmznaVLz88RWlughxtKgBUKxRll9ep3HHcAwB13LOhluQaDXm66aTavvdbMmjXbHQv01lufznEhLFkynZtvvoAlS6ZbrZlrySxcWOcI7umyXk8FSoAVCsWgyP/Z3x+7drVyyy1P8Oij1zFzprm2YTJ5cnLSSy/tY+PGw8ybN96yXGexcOEkp+2HHvo0l17a0MuatXn22b2sWbOdSy9tYPnyRpYt+5gj3iMBIaUcuJbFRRddJN9+++3T2B2FQnG6GYyA5u+3evU7bNzYwrp1f2DevPFcf/0MVqyYT2VlEJdL9Ko/Z84POXiwi8mTK3jvvS9RVubngQd+z513Pk99fTmHDplpXZcuPZ8rrpiS48stpo8f9lzONEKILVLKi/K3KwtYoRhFFCNI+X7T7H0MQ/LjH79LU1Mdv/rVbpLJNAsXTkLTMjzxxB7Wr98FQEmJj61bj7F16zGeeGIPPp+bJUvO5XOfm8E992xg3rwJBAJuDh7soro6yMGDXaxdu8OyUueyefMR1qzZzty545gypZL77ruKurqyPvvYF6dzgOxMoARYoRgF2CKaHW6VLUx2+dVXn0NLS4RPfWoq27cf51e/2sXPfrad9et38eij73D8eIz9+03RbG9PAPDrX+/miiumcPx4DwDTplVy7rnVbNhwkIqKAG63i40bW9i4sYUf/WgrBw92sW7dTlauXMyqVdeyZMl0nn12r+MWsF0LAGvWbGfFivk54gunN/RrOKFcEArFKMD+Wf/tb19OKOTliiumousGHR0JNm8+wr33vklXV5KamhBtbXFnv8bG8eze3U4slmbMmBCtrXFqaoK0tSVy2l+0qJ76+jIOHOhi9uxxbNlyhC1bjgEwZ85Y3nvvhONSmDt3HHV1ZVxySQMrVlzUpyU+UtwHpwLlglAoRgCFRKmlJcLXvvYbZs0aAwi2bz/O3//9JQQCHjo7kxw9GuXVVw8C8OSTe6itLeVb33qdj3+8jra2OHv3ttPVlcTthra2OOecU0FPT5oZM2o4fryHWCyN3+8mHPbS2grd3eZMsWDQQ3V1kEDAw9VXn0Msluab3/wkU6dW8q//utER4Ouvn8ns2eO4++5FvP76IccKf/rpvVRWBvp0EYx098GpQAmwQjHEZIuu7fuMxzVA0Ng4gS996Wl27mzN2eedd46xcOEk9uxpw+0WzJhRwyc/ORmAl18+QCSS4tVXm4lEUjQ1TSQWO04ioQMQjWqcOBHHMCTBoCkBqVSGgwfNAbF02qC01Mfq1ddx9GiUr3zleXbtamPNmu1MnlzBOedUsWLFRQQCbpJJnUDAw0MPfZqyMj8zZ47JmuorRr0L4aOiBFihOAMUsmwjkRSPPfYebW1x7rnnVXbubKWurozZs8fw2GPb2Lr1mOOLDYXc6LpE0wx8Phdut+Dw4W42bz5qtaWRTmfYu7cDAL/fjddrRiV0dSVJJHRCIQ/xuI7XayaIam01XRGXXdbA7NljaW7uprFxAjt3nmDdup20tyf44hfn4vd7WLJkek44WFmZnzvuaOLhh7f0is21yxQDowRYoThN2KK7dOn5rF79Dnfd9QI7d7YSj6epry/nt789wCuvHGTu3PFUVwd58cV9HD0apasrRVWVKdLt7QnGjg1x4oQpli4XaJrB3r0dHDtmDorV1ARzLGSfz00qlSGVylBTE2L27HFEoxpHjkQBmD17DMuWzeUPfzjBnDm1VFQEAHjwwc1ce+25fPWrTVxxxVTny8IW1kLugrNlsOx0oQbhFIqPiC20N944i3Q6w4EDXfziFzv44IMefvGLP3DOOZU0Nk7gj39sZ8+eNhIJnalTK1i8eBq7d7exefMRotF0wbZLSrwYBsTj6V7be3rMbVVVATo6klRUBPB6XbS2xpk6tYLJkyu4//6r+d3vWlix4hmWLj0fkKxbt5Obb76ANWu2O68rV17tTGAY7QNiQ4EahFMoThGalqGtLU5ra4xjx3r48Y/f4+c/38F3v/sG8+dPYNOmIzQ3d1NXV0Z1dZB9+zoZMybEjh3H0XXT4Nm/v4uf/OQ9/H63I74+n0DTzHKXCwwDR2Szqa4OMnVqJZs3H2XSpDIWLaqnszPB/fcvpra2hNtue5Y1a7bzjW9cwpw545kypRKX66Q/9oorpjouhWzXghLeM4+ygBWKLGxrdsmS6fz85zuIxVIsWFDHM8/spbzcz7p1fyAeT+N2uxg7NszUqZUcOtTN/v2dtLbGmTWrhp0723q16/O50DSj32MHgx4CATednb3z1QaDHhIJnXHjwrz44l8waVIZt9/+rJM7YdWqa/vMEKYYepQFrFDkYQvVddedx/79nfzTP73KmDEhHntsuxMTC+D1ukinDefV5vDhCLt3txGLnbRSUymdcNhDLKbniK5t5wQCbjKZDOksw9blEhiGJJHQmTNnHB/7WC0bNhxg16525swZy9SpVTz++G4AvvrVJi64YCwADz30aebPn0h+ikUV3jVyUBawYtTRlwXY3Z3koYfeJhpNcuGF47n//t+xefNRPv7xOo4ciXLoUDcTJ5aQSOh0dCT7bN92DxSDbbkWoqlpIj09Gjt2tPInf3IuPp+H++67mrq6sl5JxbPTNCqrduTRlwWsBFgxYugrlCt/mx0aNXfueL7znStIpXS2bTvOY49tc8K0GhrKaW7uZty4MJ///GyeeWYve/d25EQc9IU96OX1CoJBL9GoRn+P0fjxJVx77XSOH+9h9uxx7N7dxoIFE1mxYj6AchecBSgXhGLEUyhBi72trS3OJz85mX37Otm1q5Vw2Mu77x7jppt+iZSSZFJH08wJBmPHhpw2jx+PsWtXGxMmlLJ3bwf19eUAGIZk0qRy9u5tp6IiQCSSIhLRANOqrakJ0dBQzpYtH+T0MRh0k0hkuP76GcyZM47//M93OXSom0RC58kn/7zgeSl3wdmLEmDFsKPQMjNr1+5g0aJ6Pve5mezd284zz+xl3752Hn98DxdfPIH33jvOr3+9G4/HnKRg+2XtabVgug6iUY1oVMs53r59HUyYUArArl2mT7emJkRzczc9PeaAWySisXBhHeGwj9mzx/L9729i0aJ6tmz5gLq6MlpaIixdej733XeVk3hm7dodHDrUzcUXT+Tee688cxdQMWJQAqwYNpzM6KVz553P09OTIpXK0NzcxapVWykp8dHTY4rniy/uJx5Ps3dvBzNn1rBt2zF2724HYMaM6px2PR7QdZg5s4ZYLE1XV5KurpPC/P77nUyYUMq8eePZuvUY48eXcOxYD9/4xkJeeeUg3/3ulTzyyDvce++V1NWV8cADmwBYsKCOxYun5WT7yh4Ay56koNwLikIoH7DijNCXVbtkyXTWr99JMqmzefNR1q/fxYoVF9HaGqOlJcqmTS1UVQXp6MjNzmXmIcg4bWUzb14tBw92OgNp9fVlHDoU4TOfOZfq6iCPPvoeM2fWOK6Ho0ejrFx5tbMOWb6g2j5lO9RLhXkpBosahFMMKbaIXXzxRNavX8r69bu4887nufDCsWzbdsKpV1UVZMqUCi66qJbSUj/PP/8+O3a0OqFd48eXEImkcmaG1deX43YL6urKiMfTfOUrC3jrrSP84AebWbLkHH7wg2t54YV9OcuU2yKbL7aFUIKr+KgoAVacVrIt2mxBy95+ww3reOutI9x88wXMnz8hZ2HG0lJfjm92wYKJuFyCcePCvPnmYScm9+KLJ3LBBWN55JF3APjsZ8/l0ksb+PrXX+Tiiyfy1ltHWLXqWscHq0RTMRxQURCKj0R/VmD2irZ2boHOzgRbt35AOm3w+OO7WbnyapYsmUYmY3D33YuYONFcAaGlJYqUkkQizUMPmV/un/70NJqaJvGP//hbPv/52bS2xrn++hl4vS7uu+9qysr8XHDBOJLJNIGAlxtumMk77xxzjp/vi1UohitKgBUDkr9keL51uXr1O6xZs53GxlruvnsRl17awPPPv+/M3mpsrOWll/bz9NN7AXjhhX0Eg96cSQWRSIr6+nK2bv2A++67GoDdu9u4++5FXHbZ5F7Cf8cdC7JSIXpyVs9VFq9ipKAEWDEga9fuyLEus+Nxb7ppNhs3HgZgy5YPeP31Qyxf3khnZ4LHH9/N3LnjcmJlly6dBQhuvfVpEgmdYNDjiGZlZZB163ZyxRVTAXKWGy9EfpSBsngVIw0lwGc59lLjIFi2bG5B6zFf6LI/r127g3XrdrJ06fksXFjnlNlLlS9aVM8997zKrFk1VFQEWbZsLmBP0U3nTKwolFu2vzyzSnQVIx4pZdF/jY2NUjG6WLXqbQn3SLhHrlr1dk5Zd3dSrlr1tuzuTva5/0B17Pbz2y62fYViNAC8LQtoqrKAzzLyB9Nuumk2iUSa/PW78v2+H3Zhxf5WTFAWrOJsRwnwWUZ+PoW+1u/K9/t+WJTIKhR9owR4GHM6JgAUu4aXmkarUJx+XEPdgbOVSCTFww9v6TWNNhvbWl29+t1edfP3L6Y9OGmR2pMk+tonu55CoTg9KAt4iCiUWjEf2wq1owUSibSzcGL+/sW092H6oFAoTh9KgIeIYlwBthUaiaQIBr0kEnpO/G1/r6eqDwqF4vShckGMIFRSGIViZKJyQYwCVESBQjG6UINwCoVCMUQMygUhhGgFmk9fdxQKhWJU0iClHJO/cVACrFAoFIpTh3JBKBQKxRChBFihUCiGCCXACoVCMUQoAVYoFIohQgmwQqFQDBFKgBUKhWKIUAKsUCgUQ4QSYIVCoRgilAArFArFEPH/AW45lM0h9qOTAAAAAElFTkSuQmCC",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Plot the winner\n",
        "splot = plt.subplot(2, 1, 2)\n",
        "Y_ = clf.predict(X)\n",
        "for i, (mean, cov, color) in enumerate(zip(clf.means_, clf.covariances_, color_iter)):\n",
        "    v, w = linalg.eigh(cov)\n",
        "    if not np.any(Y_ == i):\n",
        "        continue\n",
        "    plt.scatter(X[Y_ == i, 0], X[Y_ == i, 1], 0.8, color=color)\n",
        "\n",
        "    # Plot an ellipse to show the Gaussian component\n",
        "    angle = np.arctan2(w[0][1], w[0][0])\n",
        "    angle = 180.0 * angle / np.pi  # convert to degrees\n",
        "    v = 2.0 * np.sqrt(2.0) * np.sqrt(v)\n",
        "    ell = mpl.patches.Ellipse(mean, v[0], v[1], 180.0 + angle, color=color)\n",
        "    ell.set_clip_box(splot.bbox)\n",
        "    ell.set_alpha(0.5)\n",
        "    splot.add_artist(ell)\n",
        "\n",
        "plt.xticks(())\n",
        "plt.yticks(())\n",
        "plt.title(\n",
        "    f\"Selected GMM: {best_gmm.covariance_type} model, \"\n",
        "    f\"{best_gmm.n_components} components\"\n",
        ")\n",
        "plt.subplots_adjust(hspace=0.35, bottom=0.02)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.8.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}
